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SUMMARY 


A method for the reduction of the cost of solution of large nonlinear 
structural equations was developed. Verification was made using the MARC-STRUC 
structure finite element program with test cases involving single and multiple 
degrees-of-f reedom for static geometric nonlinearities. The method developed 
was designed to exist within the envelope of accuracy and convergence charac- 
teristic of the particular finite element methodology used. 

INTRODUCTION 


At present the finite element codes in conjunction with the large, high- 
speed computers available are capable of producing reasonable solutions to 
practically all static problems conceivable in structural analysis. In addi- 
tion, well-behaved problems such as those involving small elastic deformations 
are solved relatively inexpensively and accurately. Computational difficulties 
do not arise until the stiffness of the structure becomes a function of dis- 
placement and/or displacement history. An opinion widely held is that when 
this does occur an implicit solution scheme is necessary for accuracy. All 
implicit schemes require an iterative solution where there is an attempt to 
reduce some error term to zero at each iteration. Therefore, a nonlinear 
problem is more expensive to solve and can become astronomically so depending 
upon the degree of nonlinearity and the convergence criteria used. 

In the solution of nonlinear structural equations the reformulation of 
the stiffness matrix is a first order contribution to the cost. The first 
logical step in attempting to reduce the cost would be to seek a less expen- 
sive way to update the stiffness matrix. This of course has been done with 
some success and is apparently still being researched. Looking at only the 
most recent developments or evaluations we see that Mondkar and Powell [1] 
have used the constant alpha technique to try updating the stiffness matrix 
for the modified Newton-Raphson approach. Matthies and Strang [2,3] 
have taken similar approaches bom from a paper by Dennis and More [4] on 
Quasi-Newton methods. The basic premise was that the stiffness matrix could 
be updated without going through the full process of reformulation and de- 
composition or inversion. The most popular approach was to update the stiff- 
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ness matrix by a matrix of rank two. This is known as the BFGS (Broyden- 
Fletcher-Goldf arb-Shanno) update. Crisfield [5] used a method similar to a 
BFGS update of rank one. All of these papers show conclusive evidence of cost 
reduction for certain problems. The last by Crisfield is closest in form to 
the method developed here. 

A second logical step is to reduce the number of iterations required to 
satisfy the convergence criteria. This can be done by determining an estimated 
displacement as accurately as possible. The development which follows shows 
how to do this. All forces and loads are of an incremental nature. 


DEVELOPMENT 


In general, the iterative methods of solution using the stiffness formula- 
tion will by some logical means calculate a generalized displacement for a 
given generalized load. Returning then to the elemental level, the elemental 
stiffness matrices are altered to reflect this change in shape and the total 
resistance of the structure to the applied load is determined. If the struc- 
ture is to be considered in equilibrium, the applied load must be exactly 
balanced by the resulting resistance of the structure. Any imbalance is termed 
a residual force and must be considered as an error. An attempt is made to 
reduce this error by altering the estimated displacement. The rate of conver- 
gence depends on the manner of estimated displacement selection. 

The vast majority of implicit schemes available utilize only the most 
recent residual and thereby ignore any possible trend determination. Felippa 
[6] recognized this and proposed a viable method for determining the displace- 
ment that would yield the least residual within specified limitations. This 
approach required the determination of a weighting matrix that was dependent 
upon the elements chosen and the applied loads. The development in this paper 
is independent of the physical characteristics of the elements. 

A key element in the success of the approach developed is the finite 
element method used. As mentioned before the MARC-STRUC structure program was 
used but the variational formulation of the structural equations was performed 
according to the method of Jones [7]. It is most important to have the most 
accurate method possible for the determination of the residuals. 

Considering the solution form, in Figure 1 a graph of force versus dis- 
placement is shown. The curve represents the calculated resistance of the 
structure. The original stiffness matrix, K Q , assumes linearly elastic defor- 
mation_and yields the displacement, u Q , and the residual, R 0 , for the applied 
load, F. The displacement, u Q , and residual, R Q ,_are then used in a Quasi- 
Newton fashion to update the stiffness matrix to Kj and a new displacement, u^, 
and consequently a new residual, Rj_, are calculated. Highly accurate answers 
may result, but they are clearly expensive to obtain. 

The extrapolation method presented in this paper is clearly exemplified 
by the triangle, ACE, shown in Figure 2. The method used was identical to the 
direct iteration method (shown in Figure 1 earlier) up to and through the 
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calculation of u Q , R Q and u^_, R-j_. It was the determination of the new esti- 
mated displacement, U2, that was performed differently. In_a one dimensional 
sense the residuals R Q , R^ and the distance between them, d, were used to 
calculate a scalar, a), that predicted the displacement at which equilibrium 
supposedly existed under the load, F. Of course, this was not the equilibrium 
position and a new residual, was determined. The residuals, R^, R2 > and 

the distance between them, u) d - d , were then used to predict a new equilibrium 
position. The process continued until convergence was satisfied. 

A major difficulty encountered was the determination of the scalar, a). 

In a one dimensional case it was easy enough to see that 


00 = 



( 1 ) 


However, since in general the vectors R 0 , R^ and d are heterogeneous in their 
components 1 units, a division as mentioned above is impossible even when using 
vector lengths. The solution to this difficulty came about by considering the 
units of work. In fact, this extrapolation process may be symbolically thought 
of as minimizing the work done by the residuals. In this light it was then 
decided that equating the area of the trapezoid, ABDE, plus the area of the 
triangle, BCD, to the area of the triangle, ACE, would result in an equation 
with only one unknown. Simplifying and rearranging, the following was 
obtained . 


(Ro-Ri) % d 

At this point it was decided to implement the theory and test for a one degree- 
of-freedom case and follow that with a more complex case. 


VERIFICATION 


In an attempt to determine the validity of the aforementioned extrapola- 
tion method it was determined that a one dimensional buckling problem would be 
appropriate as a first test case. The bar-spring problem of Jones was 
selected . 


Bar-Spring Problem 


In Figure 3 the dimensions used on the problem may clearly be seen. The 
length of the spring was unimportant as long as nonlinear effects did not enter 
the calculations for the spring’s deflection. The bar was modeled so as to 
allow only a change in length and no bending deformation, hence the absence of 
an El term. A load was applied at the end of the bar and spring in the direc- 
tion of deformation to render the problem one of a purely single dimensional 
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case. The buckling load was at 2.7 kg. (6 lb.) with the results tabularized 
in Table I. The exact deformation was calculated and plotted in Figure 4 to 
show the high degree of nonlinearity of the problem. 

In analyzing the results (see Table I) it was decided that a comparison of 
the values calculated against the exact values as well as a comparison of the 
number of iterations required for each method would be of use. The raised 
numbers beside the calculated displacements in Table I represent the number of 
iterations required above the original estimate to reduce the quotient of the 
calculated displacement and the estimated displacement to the tolerance indi- 
cated at the column heading. 

It should be noted that at the buckling load the tolerance required to 
obtain two significant digits accuracy, 1.001, led to a 5 vs. 26 advantage in 
iterations for the new method. However, the new method was edged by the old 
in the post-buckled region by a consistent 4 vs. 7 margin. The reason for this 
was apparently that the linear extrapolation did not follow the changing stiff- 
ness of the structure very well. If so, a better approximation would be 
obtained with a parabolic extrapolation. 

On the whole this extrapolation showed promise in this case but not of a 
clearly decisive nature. Therefore, the motivation for a more complex example 
was established. 


Ring Buckling Problem 


This problem was to determine the deflection of a ring under a uniformly 
loaded external pressure of varying values. The ring was modeled through 90 
degrees as shown in Figure 5. The 90 degree arch was broken into two substruc- 
ture. The degrees of freedom per node were 

1. Z 

2. R 

3. dZ/ds 

4. dR/ds 

with the rotations positive as shown by 0 in Figure 5. The ring was modeled 
with a modulus of elasticity of 2.1 x 10^ kg/cm2 (30 x 106 psi) and a radius 
of 51 cm (20 in.). Finally, a kicker force was applied at node 1 of substruc- 
ture 1 in the negative R direction with a magnitude of 1.5 x 10”^ kg. (3.4 x 
10“6 lb). Obviously, this was simply to force the ring into a buckled mode 
without altering the deflections due to the pressure loading. 

As there was no exact solution other than the known collapse load, the 
tolerance chosen, 1.001, was that which gave two significant digits accuracy 
for the bar-spring problem. The results obtained are shown in Table II. The 
point at which, the structure would "collapse" was 4.22 kg/cm^ (60 psi). As can 
be seen, the results were quite remarkable as the structure became softer. At 
4.18 kg/cm2 (59.5 psi) the number of iterations reached by the old method were 
not enough yet to satisfy the tolerance requirement of 1.001. The authors 
suspect that another 50 to 100 iterations would have been required. 
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CONCLUDING REMARKS 


The problem discussed in this paper was the cost of solution of large 
nonlinear structural equations. This difficulty has been and is being 
researched; however, the direction of most present research is apparently con- 
cerned with the second partial of the strain energy expression (stiffness 
matrix). This paper implies and subsequent research by the authors supports 
the supposition that the first partial of the strain energy expression (resist- 
ing force) is not being fully utilized in the determination of the new esti- 
mated displacement needed for implicit methods. It may well be determined that 
updating and/or reformulation of the stiffness matrix is occurring far too 
often in present solution techniques. 
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TABLE 1 


BAR-SPRING PROBLEM 


Load 

kg Exact 


(lb) 

Method 

(cm) 

1.1 

1.01 

1.001 

1.0001 

1.00001 

1.4 

Old 

0.59781 

0.5915} 

0.5964? 

0.5978?! 

0.5978*? 

0.59781 

(3.0 

New 


0.5915 

0.5964* 

0.5978"* 

0.5978° 

0.5978 

2.7 

Old 

2.5400 

2.0177? 

13 

2 . 4545?"'’ 

2 . 5320^ 6 

39 

2.5392^ 

2.5400?! 1 

(6.0) 

New 


2.2329 

2.5018^ 

2 . 5373"* 

2.5400 

2.5400 

4.1 

Old 

4.4821 

4.4933?" 

4.4806?" 

4.4821?- 

4.4821?- 

4.4821? 

4.4821-* 

(9.0) 

New 


4.4788 1 

4.4816 

4.4821* 

4.4821* 

5.4 

Old 

5.0800 

5.0792;" 

5.0800^ 

5.0800* 

5. 0800 * 
5.0800 

5.0800^ 

(12.0) 

New 


5.0828* 

5 . 0803* 

5.0800 

5.0800 y 

6.8 

Old 

5.4907 

5.4943?- 

5.4907;? 

5.4907* 

5.4907q 

5 - 4907 io 

5.4907 

(15.0) 

New 


5.4948 

5 . 4910 * 

5. 4907 7 

5.4907 

8.2 

Old 

5.8143 

5.8176} 

5.8146?? 

5.8143?? 

5 . 8143? 

5.8143^ 

5.8143 

(18.0) 

New 


5.8176 

5.8148* 

5.8146 7 

5.8143 


^Buckling Load 


Note: When the number of iterations is less than (2), there is NO difference between the new 

and old methods. 
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TABLE II 

RING PROBLEM (l.OOI) 


Load 

kg/cm^ 

(psi) 

Method 

Iterations 

Substructure 1 (cm) 
Node 1 
D.O.F. 2 

Substructure 2 (cm) 
Node 1 
D.O.F. 2 

.5 

Old 

2 

-.78547 E-03 

-.41397 E-03 

(7) 

New 

4 

-.78555 E-03 

-.41397 E-03 

1.5 

Old 

5 

-.25537 E-02 

-.10468 E-02 

(21) 

New 

4 

-.25545 E-02 

-.10459 E-02 

2.5 

Old 

9 

-.49439 E-02 

-.10607 E-02 

(35) 

New 

4 

-.49472 E-02 

-.10576 E-02 

3.5 

Old 

23 

-.10230 E-01 

.18172 E-02 

(49) 

New 

3 

-.10237 E-01 

.18243 E-02 

3.9 

Old 

54 

-.22451 E-01 

.12827 E-01 

(56) 

New 

3 

-.22597 E-01 

.12972 E-01 

4.2 

Old 

149* 

-.88354 E-01* 

.77978 E-01* 

(59.5) 

New 

4 

-.95669 E-01 

.85268 E-01 


^Maximum number of iterations allowed 
Convergence not yet satisfied. 




DISPLACEMENT 


Figure 1.- Direct iteration method. 



DISPLACEMENT 


Figure 2.- Direct iteration method with 
extrapolation modification. 



= 7.03 X 10« kg/cm2 (10 7 psi) 


K s =1.07 kg/cm 
(6 Ib/in.) 


-254 cm 
(100 in.) 



Figure 3.- Bar-spring problem. 


DISPLACEMENT (in.) 



LOAD 
10 (lb) 


DISPLACEMENT (cm) 


Figure 4.- Exact solution to bar-spring problem. 
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